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SUMMARY 


The task carried out under this research grant covers research on accuracy and 
efficiency of CFD strategies, error estimates for cc nvective terms, and antidiffusion. These 
basic studies (see Appendices A through C) are considered important in evaluating 
available CFD codes which will be the main activities for the next year. 


APPENDIX A 

CONVERGENCE STUDIES FOR CONVECTIVE TERMS 


Given: Burgers equation, geometry as shown in Fig. A.l. 


flu , flu 


flv 


flt 


+ u £ +v ^- 1 ' 


fl 2 U , fl 2 U 

[flx 7 fly 2 


flv . flv flv f fl 2 v , fl 2 v 

?T + U fl^ + + 


f X = c 


-fy = 0 


fx = 


1 


1 t 
1 


i 


x 2 + 2xy- rT - I 


fv = 


y 2 + 2xy- rT - f 


+ 3x 3 y 2 - 2i/y 
+ 3x 2 y 3 — 2 ux 


Exact Solution: 

u = m + x2 * 
v = j-tt + 


1 



2 


Required: 

With SUPG, Newton— Raphson iterations, and bilinear isoparametric elements, the 
coarse, intermediate and fine mesh shown in Fig. Al. Use v = 1 and v = 10 6 ; At = 10* 6 , 
1(T 4 , 10' 2 , 1, 10 2 ; T} = 0, 1/2, 1. 

Solution: 

The RMS errors vs At in log scales shown n Figs. A.l and A. 2 have trends similar 
to the linear problem (without convection terms) The error increases monotonically with 
an increase of At for t} = 1/2, whereas the error is almost independent of At for rj = 1. The 
results for = 0 are not presented as they are outside of the scale shown. If v = 1, the 
error increases as At becomes small due to round- off errors. The error decreases rapidly as 
the mesh is refined. 

Remarks 

The Newton-Raphson scheme converges after 4 or 5 iterations. It should be noted 
that the combination of the SUPG and Newton-Raphson procedures are responsible for a 
stable and accurate solution. 
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Coarse Grid 
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Intermediate Grid (Halved from the Coarse Grid) 


N N N N 


Fine Grid (Halved from the Intermediate Grid) 


N represents Neumann boundary conditions. For example, at 
node 12 for Coarse Grids » 2 xy, * x 2 , = y 2 , 

■ 2 xy are prescribed. Dirichlet boundary conditions are 
prescribed for all other boundary nodes from the exact solution. 

Use u*v*0 as the initial guess for all interior nodes. 


Fig. A. 1 Geometries for Example 6.4.1. 
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Fig. A. 2 RMS error for Example 6.4.1, v = 10^. 
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Fig. A. 3 RMS error for Example 6.4.2, v = 1. 
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APPENDIX B 

ERROR ESTIMATES FOR CONVECTIVE TERMS 

The standard Galerkin finite element method is known to have computational 
instabilities in the convection dominated flow field. To overcome difficulties involved in 
the convective nature of hyperbolic equations, the SUPG scheme may be used. The error 
estimates will be discussed for the cases with and without diffusion terms. 

Consider the convection equation of the form 

Lu = f in Q (B.la) 

u = g on T (B.lb) 

with 



The inner product of (B.la) and the test function including the sum of the trial functions 
and the numerical diffusion test functions leads tc the problem: Find u e Sh such that 
[((a • V)u,v) + h(a • V)v] - (1 + h)(u,v) 

= (f, v + h(a • V)v) - (1 + h)(g,v) V v e Sh (B.2) 

Then there exists constant c if u satisfies (B.2) such that 

H u_Q| lL 2 (ft) ^ chk+l/ 2 H u ll H r>i ( B - 3 ) 

To prove this, let u h 6 Sh be an interpolant of u satisfying the steady-state case. Denoting 
eh = u — u h and e h = u — u h , we obtain [Johnson. 1988] 

IMIJ 2 (q) = B ( e,e ) = B ( e,eh ) _ B ( e ’ §h ) = B ( e ’ eh ) 

= ((a • V)e,e h ) + h((a • V)e,e h ) + (e,e h ) 

+ h(e,(a • V)e h ) - (1 + h)(e,e h ) 

< I ll(a ' V)e|p + h-'llekll 2 + 1 ||(a V)e||* 

+ h||(a • V)e h j| 2 + \ He|p + ||e»*|p 

+ j l|e|P + h 2 ||(a • V)eh|| 2 + |e| 2 + (1 + h)|e"| 2 (B.4) 

We consider the inequality 
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Fig. B.l Typical interior nodes of (a) quadrilateral mesh and 
(b) triangular mesh for stability analysis. 
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or 


2ab < ea 2 + e^b 2 

(B.5) 

may be used in (B.4) so that 



(B.6) 

MIl 2 (A) s chk * l/ 2 NI«i...(s)) 

(B.7) 


or 


If the standard Galerkin method is used, then we have 
= ((* ' V)e*>,e h ) + (ei>,e») - (eM*) 

< ||(a ■ V)eh||2 + ||eh||2 + \ ||e»|| + ||e»|p + \ |e»| : 


H e llL 2 (!!) - ehk l! u ll H k**(n) 


(B.8) 

(B.9) 


Comparing (B.7) and (B.9), it is seen that the SUPG approximation provides 0(h 2 ) 
larger than the standard Galerkin approximation Experience has shown that, in the 
standard Galerkin methods, oscillations propagate in the crosswind and in the upwind 
directions with little damping, whereas such oscillations decay rapidly in both directions. 

The error estimates for the steady state convection— diffusion equation are evaluated 
in a manner similar to the convection equation. Consider the governing equation in the 
form 

Lu — V 2 u = f in ft (B.lOa) 

u = g on T (B.lOb) 

In view of the steady-state results and from the analysis carried out for the convection 
equation, we obtain 

||e|| L 2 (A) < ch“«/»ii«n Hk . K0) (B11) 

which is identical to (B.7). 

The mathematical justification for the gradient discontinuity test function has not 



been studied in detail. This question remaint largely empirical until further investigations 
are carried out. 
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The error estimates for the unsteady streamline diffusion method require the 
analysis of the step in (B. 7). Thus, 

nAt 

||e(t)|| < Cihk*i/2|| u || Hk+i(fi) + C2 hk J ||u|| Hk+1 ^dt (B.12) 

It should be noted that for each time step n, the analysis indicated in (B.12) is linear, but 

the solution u n will, in general, have jumps across the discrete time levels t n . 

We have discussed the error estimates as t ie convergence to the exact solution is 

being pursued. It has been shown that these error estimates are defined in various norms. 

Accuracy may be sought as high as desired by refining the computational mesh. However, 

the transient solutions are affected greatly by certain combinations of mesh sizes and 

temporal increments. It is well known that computational instability arises or amplitudes 

grow without bound if stability conditions are not met. Toward this end, we shall examine 

the so-called von Neumann stability analysis. 

Consider the convection equation 
du _ „ du 


The corresponding Taylor— Galerkin scheme can be written as 
2 4t 2 <*P1 1 Vp*' 

T'ffiTST d “ "5T 
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(B.13) 


where v n+1 = u n+1 — u n . 

For a linear element with a uniform grid, t :ie global finite element equation at the 
node j takes the form 


1 fl ,2,1 
St |5 VH + 3 Vj + F v i + i 


a 2 At 

“6~ 


Vj-i ~ 2vjj + Vj+i 
Ax 2 


fui*i 7 Uj-i] 

aJAt 

Uj-1 - 2uj + Uj+i' 

2Ax 

+ 

A** "J 


= a 
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which may be simplified as 


or 


n+l 1 . n+l n+l n+l. i 
'j + s (vj-i-2vj + vj + i) -g 

n n 

( u i t i ~ U j-l ) + 1 


aAt 

Ax 


aAt 

Ax 

aAt 


2 n+l 

( v i -i 


n+l 


n+l x 


2vj +Vj + i) 


SF 


2/ n n n 

(li j-i — 2uj + Uj*i) 


n+l 


+ 5 (1 - C2)<52 V J +1 = cA 0 Uj + \ c2<52 U j [ 1 + 5 ( 1 - c2)<52] (up - uj) 
= cAou" + 2 c2<5 2 u" 


(B. 15) 


where AoUj = (l/2)(uj+i - Uj-i), £ 2 Uj = Uj+i - 2uj + Uj-i, and c = aAt/Ax (Courant 

number). The combined spatial and temporal response of the amplitude may be written as 

u“ = ei3 x e at = e i 0j h e Otnk = e>3jhgn (B.16) 

where g = e ak (amplification factor), h = Ax, k = At, a = constant, ^ = constant 

determined by initial condition, and 
n+l n+l n ..... . 

Vj = Uj — Uj = e 1 3j^(g l)g n (B.17) 

Substituting (B.16) and (B.17) into (B.15) yields 

e i0jh(g _ l)gn + 1 (1 - C 2)(g - l)gn e i0< j-l) h(j _ 2e^ + e 2 ^) 
or 


(g - 1) [l + l (1 - c2)(e-i0h - 2 + e>0 h ) 

= ^ c(ei0 h - e-i3 h ) + ^ c 2 (e-i3h - 2 + etf h ) (B.18) 

Denoting 27 = /7h and using the trigonometric relations, it follows from (B.18) that 

(B.19) 


1 — j (1 — c 2 )sin 2 ^ 77 ic sim7 — 2c 2 sin 2 ^ j 7 


g= 1 + 

Notice that as 77 -> 0, Eq. (B.19) reduces to 

g = 1 + ic?7 - £ c 2 77 2 - 5 ic 3 77 3 + 0 (t7 4 ) 


(B.20) 
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It is clear that the stability condition requ res c < 1, which provides |g| < 1 for c = 
1. This implies that signals are propagated without distortion when the characteristics 
pass through the nodes. 

Using a similar approach, it can be shown that, for two-dimensional quadrilateral 
elements, the amplification factor at an interior node becomes (Fig. B.la) 


g = 1 — 2cxSin 2 cos 2 — 2c y sin 2 ^;jAX cos 2 ^ 
— c x c y sink x Ax sink y Ay — i c x sink x A:xcos 2 kj(Ay . 


+ CySinkyAy cos 2 


(B.21) 


where k x and k y are the wave numbers of the Fourier components, c x and c y are the Coant 
numbers in the coordinate directions, and Ax and Ay are the mesh spacings. 

The amplification factor for an interior node of a typical triangular element assumes 
the form (Fig. B.lb) 

g = 1 — 2c x sin 2 — 2CySin 2 ^(A* + c x Cy [cosk x Ax 
+ coskyAy — cos(k x Ax — kyAy)] — i j'p [sinkyAx 

A 

+ sin(k x Ax — k y Ay)] 4- [sink x Ax + 2sink y Ay 

— sin(k x Ax - k y Ay)]| (B.22) 

Now, consider a mesh with Ax = Ay and waves k x = k y , Courant numbers c x = c y . 
The stability limits for quadrilateral and triangular meshes are c < 1 and c < 1/2, 
respectively. 

Let us now examine the advection-diffusiou 

u t = -au x + i/v xx (B.23) 
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[L., 

,2 , 1 ] 

+ 3 v i + S v j+ i 

a fv] 4l - vj-i 


Vj-i - 2vj + Vj + f 

5t 

v j-i 

+ S [ 2Ax 


Ax 2 J 
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Denoting d = i/At/h 2 , Eq. (B.24) becomes 


'J + 5 ( V H ~ 2v i + v i + 0 + \ 5F “ 2 5x* ( V H “ 2v i + V J +1 ) 


Uj+l - Uj-i 


fcip ( u j-i _ 2u j + u j + 0 


or 

Vj + ^ <5 2 vj + ^ cAoVj — 2 d<5 2 vj = — cAoUj d<5 2 Uj 1 + g- A 2 — 

- ^ (-cAo + d<5 2 ) (u " 1 - u " ) = (-cAo + dA 2 )u j (B.25) 

The necessary and sufficient condition for stability according to Lax— Richtmyer [1968] 

if If - 1 6 1 - 1 (B.26) 

if |jf>0, |g| < 1 + 0(k) (B.27) 


Proceeding in a manner similar to the case of the convection equation, Eq. (B.25) reduces 
to 

e i(3jh(g _ l)gn + ^ (g - l)gn e i3< j-1) h(l _ 2(>i3*» + e 2i P h ) 

+ ^c(g - I)gn e i 0 (j-i)h( e 2 i 3 h - l) - ^ d(g - l)g n e»3< r>> h (l - 2ei3h + e 2 i 3 h ) 

= - \ cgaetf* i' 1 ’ h (e 2i P h — i) -i- dg“!ji0<j-i>h(l - 2e*3h + e 2i 0 h ) (B.28) 

or 

g n (g — 1 )e i 0 h + ^ (g — l)g n (l — 2e i & h + e 2 3h) 

+ J c(g - l)gn( e 2 i 0 h - 1) - ^ (g - 1 )g n (l - 2ei3 h + e 2 i 3h) 

= -^cg n (e 2i 3 h - 1) + dg“(l - 2eif' h + e 2i 3»») (B.29) 

Further simplifications of (B.29) result in 

g = 1 + 1+ ^ - d { cos 0h — 1) + ^ic sin/?h [-ic sin/?h + 2d(cos/?h — 1)] (B.30) 


or 
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6 - 1 = 


-ic s i n T i' TCk — 4dsin 2 


1 - 2 j-d sin 2 ^ ^ ic sin?? 


(B.31) 


Finally, as ?? -* 0, we arrive at 


g = 1 - ic?? - d + ^ c 2 ?? 2 + ic d + ^ c 2 - jlj ? 7 3 + ... (B.32) 

In comparison with (B.20), the stability criterion as given by (B.32) represents 


greater stability apparently due to the presence of physical viscosity. 
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APPENDED C 
ANTIDIFFUSION 


With all the schemes devised in the previous discussions, difficulties involved in 
steep gradients and widely disparate time and length scales may still persist. The FCT 
method originally proposed by Boris and Book [1973] and subsequently extended to the 
multidimensional case by Zalesak [1979] has been applied to the high-speed compressible 
flow problems in the context of finite elements [Eilebacher, 1984; Parrot and Christie, 
1986: Lohner, et al, 1986]. 

The FCT method is to combine the high-order scheme which may cause over- 
shooting with the lower-order scheme to stabilize through appropriate limiting processes. 
The high— order scheme may be generated from GGFE, whereas the low order scheme is 
equivalent to the SD-GGFE method which is specifically designed to stabilize the 
convective terms. Thus, the combination of FCT and SD-GGFE should enhance the 


desirability of both methods. 

The high— order scheme may follow the Standard Taylor— Galerkin process, or more 
preferably GGFE without SUPG in the present context. We write 


“ “ =-; Bi 0 (FJ,J-')-F8 j )-K o( (GJj‘-G5 j ) +HS + N a ( 

At 2 L J 

which corresponds to step 1 given by (5.3.1) , with step 2 (5.3.2) vanishing. For this 


v i i 

problem, therefore, we must modify the approach as follows: Let AU^ / At in (C.l) be 
defined between the time steps n and n+£ and write for each element, 


".US-* = \u; e + f B*F Se + HJ 


where a* is the area of each local element, 

ae = f dfl , A = f 4>< 1 8 > dfi , BJ = f *< e > d ft (C.3) 

J fie N Jfie N N J Oe N ’j 

where N represents the local node. The right-hand side of (C.2) can be evaluated locally 
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and summed through the local nodes for each element. No boundary conditions need be 
imposed at this stage. The initial guess for all variables is required to solve (C.2). 

The next step is to rewrite (C.l) in the form: 

AogUS* 1 = + 4t ( B aS F Sj :i + K^GlJj* + H"'i + NJ (C.4) 

n + 1 n + 1 

where B a p is the convection matrix with only the standard test function and Fpj 7 , G^- 2 , 

n+l 

and H a 2 are still the same as in the time step n hut will be updated as the next time step 
is incremented. The process characterized by (C. l) is seen to be equivalent to the classical 
Lax— Wendroff method. The solution of (C.4) may require several iterative cycles within 
each time increment by introducing the lumped mass matrix 

(C.5) 


(L)/ 


Aa'AUJ* 1 = EC - AS'AUj 


-A ( 9 )/ 


where EC represents the terms on the right-hand side of (C.4) except for the first term. 
The solution of (C.5) is referred to as the high— order scheme. 

The low— order scheme is devised in order I o obtain monotonic results. This can be 
achieved by adding artificial viscosity such as that, of Lapidus and Pinder [1982]. However, 
the SDM— GGFE scheme can be employed. This will require the numerical diffusion 
matrix C a( 3 and the gradient discontinuity matrix D a ^ which will be added to 
associated with Fpj in (C.4). 

A a3 AU 3 +l = At [( B a3 + c <*3 + D a3) F 3j 7 f K a3 G 3j 5 + H a 5 + N a] 
or 


= w“;l + e“4 


(C.6) 


with 


n+l n+i 

E aj 5 = ( C a3 + D a3) F 3j* 


(C.7) 


,n+i 


where E a j 2 = (C a p + DapJFpj^- The low— order s::heme is then implemented through a 
lumped mass matrix which would enhance the diffusive effect. 

aI&uJ = W a . + E . (C.8) 
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with t denoting the low-order scheme. 

Since the higher— order scheme is considered as underdiffused whereas the low-order 
scheme is overdiffused, it is necessary that we seek an optimum. To this end, subtracting 
(C.8) from (C.5) yields 

A a0 (4U§ - AU 0 ) = ( A< l s > - A$) AUj - E; aj (C.9) 

with h denoting the high-order scheme (C.5). He e, it is seen that an appropriate limiting 
or antidiffuse process for (AU|j - AUg) would be crucial To prevent undershoots or 
overshoots. The combined high and low— order schemes may be written as 

US ' 1 = us + u£ + (AUS - au£) 


or 

AOS* 1 = U* + >)AU 0 


(C.10) 


where T] is the limiting parameter with the ranges 
— 1 < rj < 1 


Ua =US+U* 
AU„ = AUS - auJ; 


(C.ll) 

( 0 . 12 ) 

(C.13) 


Here, the limiting process must be carried out at the element level in order to ensure local 
conservation requirements by means of the limiting parameter rj. The magnitude and sign 
of the limiting parameter q can be determined as iollows: 

(1) To determine the nodal values in the form of 

US*' = u‘ + e l(?«AU e ) (C.14) 

(2) First determine rj e 

, e = mmf R * if ’« 4U « >0 (C.15) 

R' if J?eAU e < 0 

where 


4 


R- = 


min (1,Q‘ /p") 
0 


if 

if 


p + > 0, p' <.: 0 

P ! = 0 


(C.16) 
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where p' is the sum of all positive (negative) antidiffuse element contributions to node o, 

1 n | max /a \ 

P = 1 min ) (0 ’ ”' ) 


(C.17) 


and Q' is the maximum (minimum) increment (decrement) node a allowed to achieve 
in (C.4), 


O' = u max - u 

a a 

m in 

where U£ ax can be determined as follows: 
(1) Maximum (minimum) nodal values U£ 

T e 


(C. 18) 




max 

min 


(C.19) 

(C.20) 


max 

Tjm in _ 


max 

min 


(K ■ K) 

(2) Maximum (minimum) nodal values of element 

U« = { S } < u .' «»-• V 

where 1,2,3...N represent the nodes of element e 

(3) Maximum (minimum) U e of all elements surrounding node a 

(U t , U 2) U 3 , ... U M ) (C.21) 

where 1,2,3, ... M represent the elements surround ng node a. 

With (C.21) substituted into (C.18) and subsequently to (C.16), (C.15), (C.14), and 
finally to (C.10), we complete the limiting process 

In summary, the limiting process should generate no new maxima and minima in 
the solution, nor should it accentuate already existing extrema. Such a prescription 
obviously maintains positivity. To this end, we must correct the antidiffuse mass fluxes. 
Note that the antidiffuse fluxes are limited term by term so that antidiffuse flux transfer 
can push the flux value at any node beyond the flux value at neighboring nodes. This is 
the origin of the name Flux Corrected Transport [Boris and Book, 1971]. 

The limiting process described above may be applied to a single variable. Density 
is the most logical choice in compressible flows. Although, this will reduce the amount of 
computing time, adequacy of involving other variables also should be verified in each 



problem under study. Furthermore, the most acceptable limiting procedure for antidif- 
fussion remains an open question, subject to extensi ve future research. 



